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Abstract 

The temperature dependence of the topological susceptibility in QCD, Xt, essentially 
determines the abundance of the QCD axion in the Universe, and is commonly esti¬ 
mated, based on the instanton picture, to be a certain negative power of temperature. 
While lattice QCD should be able to check this behavior in principle, the temperature 
range where lattice QCD works is rather limited in practice, because the topological 
charge is apt to freezes at high temperatures. In this work, two exploratory studies are 
presented. In the first part, we try to specify the temperature range in the quenched 
approximation. Since our purpose here is to estimate the range expected in unquenched 
QCD through quenched simulations, hybrid Monte Carlo (HMC) algorithm is employed 
instead of heatbath algorithm. We obtain an indication that unquenched calculations of 
\t encounter the serious problem of autocorrelation already at T ~ 2 T c or even below 
with the plain HMC. In the second part, we revisit the axion abundance. The absolute 
value and the temperature dependence of Xt in real QCD can be significantly different 
from that in the quenched approximation, and is not well established above the critical 
temperature. Motivated by this fact and precedent arguments which disagree with the 
conventional instanton picture, we estimate the axion abundance in an extreme case where 
Xt decreases much faster than the conventional power-like behavior. We find a significant 
enhancement of the axion abundance in such a case. 



1 Introduction 


It is widely believed that the instanton calculus in QCD [1] makes sense at high temperatures. 
The asymptotic freedom ensures the perturbative expansion sensible and, most importantly, 
the infrared divergences from the large instanton contributions are cut-off by the Debye 
length, providing finite results after the integration over the instanton size (2j. 

In the semi-classical instanton picture, physics becomes ^-parameter dependent by the 
instanton contributions to the path integral [ 3 ]. The instanton calculus indicates that such a 
dependence is proportional to the product of the quark masses, rri^ f and Aq CD , where b is the 
beta-function coefficient, b = llN c /3 — 2Nf/3. For example, the topological susceptibility, 
Xt = (<9 2 / d9 2 )V e ff(0) , is proportional to m^ f AQ CD T 4 ^ N f~ b by the dimensional analysis. 

In the QCD axion model to solve the strong CP problem, the 0 angle is promoted to 
the axion field a(x)/f a , where f a is the axion decay constant [31 E, Q Eli IE, 2 HE, II] • The 
topological susceptibility is directly related to the mass of the axion as Xt = m afl- Therefore, 
the temperature dependence of xt discussed above represents that of the axion mass, which 
is important for the calculation of the axion abundance in the Universe. In the misalignment 
mechanism for the axion generation in the early Universe, the axion number density is 
proportional to the axion mass at the temperature at which the axion field starts coherent 
oscillations mmm- The instanton based estimation of the temperature dependence is 
commonly used in the literature, and predicts that the axion can naturally be dark matter of 
the Universe when m a ~ 10 -5 eV, whereas the astrophysical bound on m a is m a < 10 -2 eV. 
(See, e.g., [IB]-) The allowed region, 1CT 5 eV < m a < 1CU 2 eV, is called the “axion window.” 

There have been arguments which imply the disappearance of the instanton effects at 
high temperatures. In Ref. |T6j, it has been argued in two-flavor QCD that the difference 
between the two-point functions of iso-singlet and iso-vector scalar operators vanishes in the 
chiral limit as fast as 0(?n 2 ) when the temperature is higher than the critical temperature. 
By the Ward-Takahashi identities, this immediately means that xt should be at most an 
0(rriq ) quantity, that is different from 0(m 2 ) from the instanton calculus. Moreover, in a 
recent paper m, it is claimed that xt is 0(m^) with an arbitrary N and thus it is vanishing 
even with finite quark masses. Although these are somewhat surprising results, it is certainly 
possible that the instanton picture fails to describe the full quantum mechanical vacuum [T8]. 
If that is the case, the estimation of the axion abundance is significantly affected. 

The lattice simulations can in principle discriminate whether instantons make sense or 
not. However, the current situation is not conclusive. In Refs. H2 122, the temperature 
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dependence of various susceptibilities are measured, which seem to be consistent with the 
semi-classical instanton picture and but at the same time not to exclude other possibilities. 
Meanwhile, the analysis in Ref. (21j suggests the effective restoration of Ua{ 1) symmetry 
right above the critical temperature, indicating the disappearance of the instanton effects. 
The presence or absence of the Ua( 1) symmetry at the critical temperature T c can also 
be inferred through the nature of the chiral phase transition of massless two flavor QCD as 
attempted in, e.g., Ref. [22] (For a related work analyzing the system with the renormalization 
group flow, see Ref. [23]). See Ref. pHj for a new proposal to extract the instanton effects in 
lattice QCD, and Ref. [25] for the approach from holography. 

One unavoidable problem on the lattice is that the range of the temperature in which 
one can reliably study the temperature dependence of xt is rather limited, because the net 
susceptibility XtV, where V is the volume, rapidly decreases with temperature and hence the 
topology tends to freeze at high temperature. This problem is fatal especially in dynamical 
lattice QCD simulations since one can not realize arbitrary large volumes while keeping both 
the light quark masses and lattice artifacts reasonably small. Then, it is natural to ask to 
what temperature xt can be reliably calculated in dynamical QCD with a typical setup. From 
this viewpoint, determining the precise value of quenched Xt{T) with large volumes and huge 
statistics would not provide useful information. 

For the reliable calculation of Xt- it is preferable to use lattice chiral fermion for dynamical 
quarks. But such dynamical simulations at high temperatures are too costly to start with. As 
the second best, we choose to perform quenched simulations with hybrid Monte Carlo (HMC) 
algorithm and Iwasaki gauge action. The reasons for these choices are as follows. To answer 
the above question, the simulation environment needs to be as close to that in dynamical one 
as possible. Long auto-correlation of the topological charge is one of the serious bottlenecks in 
dynamical simulations, and the heatbath algorithm, which is available only in the quenched 
approximation, is faster but has totally different property from HMC in this regard. Iwasaki 
gauge action is known to result in relatively long autocorrelation time for the topological 
charge [251 [27]. Thus, we expect that more information about dynamical simulations will 
be gained from quenched one by taking HMC and Iwasaki gauge action rather than taking 
heatbath and the standard plaquette gauge action. 

Another issue to be addressed is the definition of topological charge Q. There is a subtlety 
in measuring Q if one adopts the field theoretical (or bosonic) definition. This method 
requires a suitable number of coolings of gauge configurations. If one applies it too much or 
too less, one would miss the right value of Q. Even if one ceases the cooling adequately, Q 
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thus obtained is not an integer, and in the worst case it may be right in the middle of two 
integers, which can cause misidentification of Q. Misidentifications are potentially dangerous, 
especially at high temperatures ( i.e. when ( Q 2 ) = Xt(T) V is tiny), because it can significantly 
affect Xt(T) even if it occurs only rarely. Importantly, it is not possible to know the right value 
of Q without comparing that obtained with the fermionic definition based on Atiyah-Singer 
index theorem. From this viewpoint, the use of the fermionic definition has an advantage, 
although it is much more demanding than the bosonic one. Studying the autocorrelation 
of Q usually requires high statistics. Thus, this choice of the definition seriously limits the 
lattice size to small. 

In this paper, we first explore to which temperature quenched simulations using HMC and 
Iwasaki gauge action are able to obtain reliable results for xt , where the topological charge Q 
in each configuration is determined by the number of zero modes of overlap Dirac operator. 
Because of this choice, the lattice volumes are limited to relatively small. Our quenched 
calculations show that xt is undetermined above T = 2 T c for Nt = 4 and T = 1.5 T c for 
Nt = 6, where Nt is the number of lattice sites in the time direction. 

This observation suggests two possibilities. One is that the result with Nt = 6 is right 
and xt suddenly decreases at T < 1.5 T c . In the language of the semi-classical instanton 
picture, this behavior can happen if there is a sharp short-distance cut-off in the instanton 
size parameter p so that small instantons do not contribute to the path integral. (See [28) [29] 
for the instanton based models and lattice results suggesting this picture.) If xt exhibits such 
a sharp fall off, the estimation of the abundance of the QCD axion in the Universe may be 
significantly affected. However, this possibility appears to be unlikely because our value of xt 
with Nt = 4 is reasonably consistent with the previous calculations to T ~ 1.75 T c . Another 
possibility is that, in the standard HMC with a fixed acceptance ratio, the suppression of 
Xt due to the long autocorrelation, which increases with a lattice volume, overcomes the 
enhancement by enlarging the volume. Through the analysis of the autocorrelation, we 
confirm the latter to be the case. 

Outcome of the first part is that with HMC the reliable calculation of xt is difficult for 
T > 2 T c even in the quenched approximation. Since quenched simulations are always easier 
than dynamical one, the above outcome or even worse should hold for dynamical QCD. 

The second part of the paper goes as follows. We first emphasize that even in the quenched 
theory it is hard to justify instanton calculus from the first principles because instanton 
description consists of the perturbative expansion in terms of the strong coupling constant 
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a s (T ) and its prediction is only reliable above T = a few GeV. The same is true or the 
situation is even worse for dynamical QCD. Furthermore, there exist convincing arguments 
that the presence of light dynamical quarks could drastically change xt above T c as stated 
above. Therefore, not much has been known about the unquenched theory, and there are 
many degrees of freedom in the choice of T dependence of xt above T c . 

We revisit the axion abundance with some extreme but yet allowed temperature depen¬ 
dence of xt- We find that the axion abundance is significantly enhanced compared to the 
conventional estimates with non power-like behavior of xt and, for some cases, the axion 
window is closed depending on how quickly Xt disappears. 


2 Topological susceptibility on the lattice 


In the continuum, the topological susceptibility, Xt, is defined as 

2 


Xt = (jTj) J i i x^F^{xy*'^°'FX.FXm- 


This quantity is related to the axion mass, m a , through 

2 j-2 

Xt = mj a , 


( 1 ) 

(2) 


where f a is the axion decay constant. Through the Atiyah-Singer index theorem, one can 
also find 


Xt = 


(Q 2 ) 

v 


(3) 


where the topological charge, Q = n + — n_, is the difference between the numbers of the left 
and the right-handed zero modes in the eigenvectors of the Dirac operator. 


On the lattice, the definition of gauge field strength tensor, F ^ w , is not unique, and so is 
Xt if one simply follows Eq. ([!]). With the definitions of commonly used in the literature, 
Q measured by 

Q = (d?) J d ‘ x ^ F ^’ w 

takes non-integer values in general. Thus, a technique called “cooling” is usually applied to 
guess the right integer for Q. 

Similarly to the continuum, an alternative is to count the zero eigenvalues of the lattice 
Dirac operator. This requires lattice Dirac operators to satisfy the Ginsparg-Wilson rela¬ 
tion m- The overlap Dirac operator, D ov , is known as such an operator m, with which 
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one can express Q as 


Q = TiT 5 , (5) 

where 

r 5 = 75 ~ = \[lb- sgn (H w (-M 0 ))} . ( 6 ) 

Here, H\y = 75 -Du/(—Mo), and Dw(-Mq) is the Wilson Dirac operator with a negative 
mass parameter, —Mq. This definition provides an integer value of Q in each configuration. 
The trace of the spinor indices of Ts, trT 5 , can provide a definition of the local value of the 
topological density, (1 /64n 2 )e tJ ’ l/pa F^F^. 

As stated in the introduction, the bosonic definition has the chance to misidentify Q. 
While, below and around the (pseudo-)critical temperature, previous works have revealed 
that the bosonic and fermionic definitions of Q give consistent results for Xt, it is n °t well 
established above 1.5 T c yet. Especially, we have to keep in mind that a tiny amount of 
misidentihcation of Q may bring significant effects to the resulting xt at high temperatures 
as follows. 

Whether the instanton picture are correct, it is certainly true that xt decreases with the 
temperature and xtV as well. Since xtV represents a net width of the fluctuation of Q , 
the fraction of configurations with non-zero Q becomes much smaller than that with the 
vanishing Q when XtV 1- Suppose that lV con f configurations are generated on the lattice 
volume of Ny = Nf x Nt and that q of them belong to either Q = +1 or Q = — 1 while the 
others to Q = 0. This yields a^xt = {Q 2 )/N v = (q + S mis )/(N conf N v ), where 5 m is/N con{ 
represents the rate of misidentihcation. Then, it should be noted that, even if 5 m i s /N con f 
is tiny, arxt may significantly deviate from a real value because the fraction of nonzero Q 
configurations, q/N c on f, is also tiny at high temperatures. To avoid the misidentihcation, we 
adopt the fermionic definition in Eq. © for the evaluation of xt- 

3 Xt(T) in the quenched approximation with HMC 

3.1 lattice parameters 

To calculate the topological susceptibility at finite temperatures, we perform lattice simula¬ 
tions of 517(3) pure Yang-Mills theory around and above T c . We employ the Hybrid Monte 
Carlo (HMC) algorithm to study the autocorrelation of topological charge in HMC for the 
reason described in introduction. To avoid the ordinary hnite size effects, the aspect ratio is 
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V/a 4 

P 

T/T c 

H of traj. 

5t 

acc. 

Xt/Tc 

16 3 x 4 

2.288 

1.00 

5000 

1/55 

0.79 

0.126(12) 


2.450 

1.34 

15000 

1/60 

0.78 

0.176(17) xlO" 1 


2.522 

1.50 

19970 

1/65 

0.78 

0.66(9) xl0~ 2 


2.623 

1.75 

39860 

1/62 

0.77 

0.43(11) x10“ 2 


2.716 

2.00 

33660 

1/67 

0.79 

0.90(41) x10 -3 


2.802 

2.25 

43250 

1/65 

0.79 

[0.9(9) xlO -4 ] 

18 3 x 6 

2.445 

0.89 

4500 

1/80 

0.82 

0.148(18) 


2.515 

1.00 

19500 

1/80 

0.82 

0.143(14) 


2.710 

1.34 

17560 

1/80 

0.80 

[0.57(21) xlO" 2 ] 


2.794 

1.50 

15710 

1/80 

0.79 



3.018 

2.00 

17050 

1/80 

0.78 


24 3 x 6 

2.445 

0.89 

4600 

1/77 

0.73 

0.111(15) 


2.794 

1.50 

2140 

1/100 

0.80 

[0.28(19) xlO” 2 ] 


3.018 

2.00 

4030 

1/100 

0.78 



Table 1: Simulation parameters. The trajectory length is always set to one. acc. denotes the 
acceptance ratio. Topological charge Q is calculated every ten trajectories. The statistical 
errors are estimated by the jackknife method with the bin size of 500 trajectories. The blank 
indicates that xt can not be calculated due to too rare Q changes. 


set to three or four, which is usually considered to be safe. In addition, we have to recall that 
the size of fluctuation of Q is directly affected by a volume through (Q 2 ) = xt V, which is a 
finite size effect specific to and important in the calculation of xt- The calculation is carried 
out on three lattice volumes using the Iwasaki gauge action, 16 3 x 4, 18 3 x 6 and 24 3 x 6. 
The topological charge Q is calculated at every ten trajectories. The statistical errors are 
estimated by the standard jack-knife method with the bin size of 500 trajectories. 

The correspondence between (5 and T/T c for the Iwasaki gauge action is read from 
Ref. [32]. We accumulated 2000 to 40000 trajectories, depending on the simulation pa¬ 
rameters, which are summarized in Tab. [l] 

Lattice volumes chosen in the present work are somewhat smaller than other pure YM 
calculations. The reason is as follows. In order to reduce the ambiguity associated with 
the definition of Q as much as possible, we decided to use the fermionic definition for Q 
to calculate xt- The purpose of this work is to find the highest temperature to which the 
reliable extraction of xt is feasible in HMC, or in other words, at which temperature the 
autocorrelation time becomes unreasonably long. Such a study usually requires high statistics. 
Furthermore the fermionic definition is much more expensive than the bosonic one. Thus this 
choice of definition constrains lattice volumes to be small. 

Here let us mention the difference of this work from Ref. [33], which studies the T 
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dependence of xt around the critical temperature, using the fermionic definition for Q. Two 
major differences are in the range of the temperature explored and rigorous application of 
the overlap Dirac operator. The former is originating from the different motivation. Let 
us explain more about the latter. In Ref. (33] . the identification of Q is made using an 
approximate solution of the Ginsparg-Wilson equation, while we stick to the exact one. They 
differ in important way. To check the approximated solutions, in Ref. [33] the chirality of the 
zeromodes (V’JtsV’i) is examined, and found that its absolute value ranges from 0.4 to 0.9. 
With our exact method, it only takes +1 or -1 within a rounding error. Therefore, in Ref. (33] . 
the possibility of the misidentification cannot be excluded thoroughly. Indeed, they observe 
that 2 to 9 % of configurations are misidentified by comparing their approximated solutions 
and the exact one below the critical temperature. Importantly, it is not clear how well the 
method based on the approximation works at higher temperatures. In our calculation, the 
size of zeromodes and near-zeromodes of Hermitian overlap operator H ov differ by a factor 
of O(10 6 ), and we checked for all zeromodes whether the complex conjugate pair is absent. 
Thus, no ambiguity exists. While 2 to 9 % of misidentifications will not significantly affect 
Xt at low temperatures, it does at high temperatures as x.t itself is very small there. Note 
that both of the two differences are essential in the study of axion. 

3.2 autocorrelation of Q 

Figure [I] shows the history of Q at several simulations. It is seen that, on our smallest lattices, 
Q changes reasonably frequently below T < 1.75 T c . But at T = 2 T c , we only observe Q = — 1 
for nonzero Q. Further increasing T to T = 2.25 T c , non-zero Q was only observed in one 
configuration in spite of our largest statistics. Thus, xt is less reliable at T > 2 T c . Another 
remark, which is not plotted, is that while \Q\ fluctuates over the range of —8 to +8 and —2 to 
+2 at T = T c and 1.34 T c , respectively, Q only takes a value of —1, 0 or +1 at T = 1.5 T c and 
1.75 T c . While this is expected as xt decreases, the central values and the statistical errors 
for those temperatures needs to be checked by comparing other works. (The comparison is 
made later in Fig. [3j) 

It is naively expected that Q fluctuates more frequently at larger lattices. However, the 
plots in Fig. [l] shows the opposite tendency, which indicates that the suppression of xt due to 
the long autocorrelation in HMC with a fixed acceptance ratio overcomes the enhancement 
by enlarging the volume for our choice of lattice setup. 

To quantify the autocorrelation, we estimate the integrated autocorrelation time, r; nt . 
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Figure 1: History of topological charge Q above T c . 

Figure [2] shows that Tint takes 10 to 20 trajectories at T < 1.5 T c on 16 3 x 4 and increases 
to about 40 trajectories at T < 1.75 T c . Note that we could not obtain r^t at T > 2.0 T c . 
We observe the qualitatively similar behavior for 18 3 x 6 lattices except for the T = T c data, 
which shows T m t larger than other results obtained at T ~ T c . 

3.3 T dependence of \t 

The figure [3] shows the T dependence of x.t- It is seen that, restricting the data points to the 
one in which Q shows a recognizable fluctuation, the results obtained from such fluctuations 
are consistent with the existing results, for example, given in Ref. |33l 1.341135j . It is important 
to note that we have only observed \Q\ < 1 at T = 1.50 T c and 1.75 T c , but the resulting yys 
nevertheless reasonably agree with the previous results. It is also found that xt with different 
volumes shows consistency within two standard deviations as long as the data showing a 
recognizable fluctuation of Q are concerned. 
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Figure 2: Integrated autocorrelation time for 
lattices (right). 



Q on 16 3 x 4 (left) and on the other larger 


To see the consistency of our results with the instanton calculus, we fit the five results in 
the range of 1 < T/T c < 2 to the form of oc 1 /T 7 , which is suggested in the leading order 
instanton calculus. The fit quality is reasonable (x 2 /d.o.f.= 0.97). 


From these observations, we realized that with HMC it is difficult to obtain the reliable xt 
above T = 2 T c or even 1.5 T c depending on the lattice volume, which immediately indicates 
the difficulty in dynamical simulations. Even if much faster computers were used, this upper 
bound will not change significantly. Thus, to estimate the T dependence of xt at 0(10T c ), 
we have to make a long extrapolation using those obtained in such a rather limited range of 
T. To push the limit upward as high as possible, it is crucial to explore the HMC parameters 
or improve the algorithm. 


Inspired by Refs. [MISIj, we tried, as an attempt, to enhance the number of configurations 
with nonzero Q by inserting 


X = det 


( K + v 2 

V + e2 



(7) 


to the path integral, where is a positive integer. Then, xt is calculated through 


<(Q 2 / V)XX~ 1 ) _ <(Q 2 / V)X~ 1 )x 
{XX-1) (X-I)x 


( 8 ) 


where (■ ■ - )x denotes the average over the configurations generated with the extra reweighting 
factor X. For ^ > e, the insertion of X enhances the eigenvalue density in the small eigenvalue 
region, whereas eigenmodes with eigenvalues A n are left untouched. Since, when the 
topology changes, the smallest eigenvalue of H\y passes through zero, the above factor is 
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Figure 3: The temperature dependence of the topological susceptibility, xt in unit of T f. The 
results of 16 3 x 4 (filled circles), 18 3 x 6 (filled squares) and 24 3 x 6 (filled triangles) lattices 
with the standard HMC are shown. The dotted curve is the fit of five data points of 16 3 X 4 
lattices obtained in the range 1 < T/T c < 2. 


expected to increase such opportunities. However, after we performed some trial calculations, 
we realized that this method does not always work and the fine tuning of /it, e and are 
required. Further investigations to improve the situation is in progress. 


4 Effects of dynamical quarks 


Let us discuss what would happen when we include the dynamical quarks. The naive guess 

Nt . Nt 

would be that xt hi the Yang-Milles theory is multiplied by a factor of m q 1 /Aq^ d since xt 
should vanish when one of the quark masses goes to zero. 


There can be more drastic possibilities. If we accept the claims of the axial 17(1) restora¬ 
tion in two-flavor QCD dsnm, the 0(m q ) contributions to xt is forbidden in two-flavor QCD. 
Therefore, the possibility of just multiplying by m q f / Aq^ d is not consistent. The results 
of Ref. |17| even forbid contributions with any power of rn q for a small m q . An extreme 
possibility one can consider is 


XtiT) 


m Q^QCDi 

m 9 A QCD e 


- 2c{m q )T°-/T* 


T < T c , 
T > T c , 


(9) 


with c(m q ) —> oo as m q —> 0, so that xt cannot be expanded around m q = 0. Note that 
the results of Ref. m is contained as a special case of eq. ([9]). Since no unquenched result 
of xt is available at high temperatures, we take c(m q ) as a free parameter in the following 
discussion. 
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5 Axion abundance 


We discuss the impact on the axion abundance in the Universe for the case where the dilute 
instanton gas approximation fails at high temperatures. One of the source of the axion energy 
density today is the coherent oscillations of the axion field started in the early Universe. The 
equation of motion for the axion field, o, is given by 

a + 3Ha =—m 2 (T)a, (10) 


where H is the Hubble parameter and m a {T) is the axion mass at a temperature T. The 
axion mass is related to xt as in Eq. <©• This equation of motion leads to an oscillating 
solution. Although the axion mass is temperature dependent, and thus time dependent, the 
axion number density divided by T 3 , n a /T 3 , stays constant during the oscillation and its 
value is given by 


T 3 


m a(T*)f 2 9 2 

T 3 


( 11 ) 


where T* is the temperature at which the axion field starts to oscillate, and 9 is the initial 
amplitude of a/f a . In the conventional estimates, one uses m a {T ) oc T~ 4 derived by the 
instanton gas approximation. The temperature T* is obtained by equating m a (T ) and the 
Hubble parameter, H(T ) ~ T 2 /Mp \, with Mp\ being the Planck scale, such as 


1 GeV 


m„ 


1/6 


( 12 ) 


. 10 -5 eV/ 

By using this value and setting m a (T*) ~ 3id(T*), the energy fraction today is finally given 
by 


n ~ o 2 • e 2 ■ ( m ° \ 

“ y 1^10-5 eV/ 


-7/6 


.10-eW . (13) 

for 6 <C 1, and it is about a factor of two larger for 9 ~ 0(1) [38( [39] • According to this 
formula, the axion can naturally be the dark matter of the Universe for m a ~ 10 -5 eV which 
corresponds to f a = 6 x 10 11 GeV. (See Ref. [30] for a recent detailed calculation.) 


If the axion mass m a {T) decreases much faster as discussed in the previous section, 
the axion starts to oscillate near the critical temperature T c ~ 150 MeV independent of 
the axion mass. More importantly, the axion number density gets fixed when the time 
scale of the oscillation, l/m a (T), is comparable to that of the change of the axion mass, 
(rh a (T) / m a {T ))~ l , so that the change of the mass is adiabatic, i.e., 


m a {T) 


1 dmn a (T) 


din m a (T) 

m a (T) dt 

n ) 

dlnT 


(14) 
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The condition is automatically satisfied when m a (T ) ~ H(T ) and m a {T) oc T n with n ~ 0(1), 
but not necessarily true for exponential functions. 

In the model discussed in the previous section, m a (T) is given by 

m a {T) = m a {T c )e- c{T ~ J - T ^ /T ° , (15) 

where the c parameter may be much larger than 0(1). The oscillation starts when both sides 


of Eq. (14) become comparable, and thus 


m a (T*)~2cH(T*)(^j . (16) 

Since m a (T) is a steeply varying function, we expect T* ~ T\, where T\ is defined by 

m a (Ti) = 3 H{Ti). (17) 

This causes an enhancement of the axion density by a factor of about c(Ti/T c ) 2 compared to 
the conventional estimate in addition to the enhancement due to the shift of the oscillation 
temperature. 


A numerical solution of the equation of motion in Eq. (10) shows that the axion number 
density is approximately given by 


n a (T ) _ 0.5c x 377 (Ti)/ 2 0 2 T 2 _ 0.5c x 3 H{T c )f 2 9 2 T x 

' T 2 


T 3 


Tf 


T 3 


(18) 


T\ 


1 + * [ 12 + log 


(19) 


Here 

m a (T c ) \ Y ] 1/2 
10 _3 eV ) ) 

Here we used 3 H(T C ) ~ 5 X 10~ 20 GeV, and assumed (Ti/T c ) 2 <C m a (T c )/3H(T c ). Putting 
it all together, the axion energy density today is given by 

1/2 


Si„~O.2-0 2 - (jAAy) 1 X 2 - 5c 


1 + ^ ( 12 + log 


m a (T c ) 
10 -5 eV 


( 20 ) 


Even for c = 0(1) and m a = 10 -5 eV, we find that the axion energy density is enhanced by 
an order of magnitude compared to the conventional estimate. The enhancement is larger 
for larger c. The axion window is closed when c ~ 400. For an extremely large value of 
c, i.e., c > 4 x 10 5 • (m a / 10 -5 eV), the axion starts to oscillate at T = T c . In that case, 
G a ~ 2 x 10 5 • 8 2 independent of c. 

We have also examined the case of xt oc exp(—2 cT/T c ) rather than xt oc exp(—2cT 2 /T 2 ) 
in Eq. ©■ We hnd a similar enhancement factor 1.0 c instead of 2.5 c in Eq. (|20[). 
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In the scenario where the Peccei-Quinn symmetry breaking takes place after the inflation, 
one should take an average value of 9 , (O' 2 ) ~ n 2 . Also, in addition to the misalignment 
mechanism considered above, there are contributions to the axion density from the axionic 
strings and domain walls mmmm- Both contributions receive an enhancement due 
to the delay of time scale of the axion production which leads a milder dilution due to 
the cosmic expansion. The enhancement factor of fl a compared to the estimates based on 
m a (T) oc T -4 , Q c a onv -, is the ratio of the temperature in Eq. (12) and the temperature T\ at 
which the adiabatic condition is satisfied. The enhancement factor is 




~ 6.6 x 


m n 


10 -5 eV 


1/6 


1 + * ( 12 + log 


m a (T c ) \ 
10- 5 eV ) 


- 1/2 


( 21 ) 


string+wall 

For m a = 10” 5 eV, for example, the enhancement is a factor of a few for c = 0(1). The value 
I string+wall is found in Ref. @5]: 

-7/6 


of n c „ 


n c „ 


4 x 


string+wall 


m Q 


10 -5 eV 


( 22 ) 


6 Summary 

The temperature dependence of the topological susceptibility xt essentially determines the 
abundance of the QCD axion in the Universe, and can be calculated using lattice QCD in 
model independent way. However, at high temperatures, the lattice calculation suffers from 
various difficulties such as the numerically demanding definition, possible misidentifications 
and the long autocorrelation of the topological charge. 

As the first step towards the precise calculation, we have performed the lattice calculation 
of xt in the SU(3) Yang-Milles theory, focusing on the autocorrelation time in HMC, and 
found that the problem starts to show up around T ~ 2.0 T c , This immediately indicates 
the difficulty in dynamical simulations at the similar temperature while estimating the axion 
abundance in the conventional scenario requires the T dependence up to 0(10 T c ). To fill 
this gap, we need exploratory studies of HMC parameters or even algorithm itself. 

Meanwhile, in order to see the importance of understanding the temperature dependence 
of xt for the axion abundance in the Universe, we have calculated the axion energy density, 
f2 a , in an extreme case where xt drops exponentially above the critical temperature. It 
is found that the non-adiabatic growth of the axion potential makes the axion abundance 
significantly larger, possibly closing the axion window. In this particular case, the behavior 
of xt around T c becomes important. 
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Our findings provides strong motivation to further serious studies of topology on the 
lattice and instanton at medium temperatures. 
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Note added 

While this paper was being completed, Ref. [56] appeared which estimates xt at high tem¬ 
peratures in the quenched approximation with the heatbath algorithm on large volumes and 
uses it to discuss the axion abundance. It turns out that our results at Nt = 4 are consistent 
with theirs to T ~ 2 T c , which give us a more confidence that the autocorrelation time of 
topological charge increases with volume in HMC and hence increasing volume does not 
resolve the problem. 
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